「R | parallel」第一篇｜gene level的FST计算

2023-05-15, 700 words, 4 min read

写在前面

聊聊我自己。群体遗传学研究的对象是allele，而针对一个population，就算是一个population，其内部存在的variants也是相当多的（由于genetic drift的影响），所以要处理的下游数据难免就会呈现records多到“上天”的情况

‍

R | parallel的使用

虽然R语言是一门针对统计学的语言，但并不就代表它不可以实现最基本的计算功能，multiprocessing。

需要使用的packages，parallel。

就拿我自己做的分析来举例：

现存一个vcftools计算得到的FST结果文件，分别记录了从chromosome1到chromosome21的per-site FST计算结果；
我的reference gene BED从UCSC Genome Table Browser下载得到，从共记录了两万多条gene records

如何计算得到每一个gene的FST？

前提条件说了很多了，只是为了希望做分析的人能够理解这个问题。

但是在给出自己的代码之前，先给出R parallel的基本使用。

R | parallele | multiprocessing

直接使用multiprocessing方法时，不需要自行创建cluster来预设对应包含线程数的sockets，而用起来就像lapply一样，但是与其不同的是，lapply是将所有的函数依次应用到输入的每一个list上，而multiprocessing的lapply则是同时发动所有的任务，

使用函数，

mclapply()
mcmapply()

R | parallele | sockets

sockets类型的multiprocessing与上述提到的基本方法，在使用上的不同需要首先创建一个clusters，同时在完成运算之后需要释放这个cluster，
同时还需要将应用的数据export到cluster中

cluster_1 <- makeCluster(16)  # 可以使用detectCores()来检测总共有多少个核
stopCluster(cluster_1)

‍

R | 我的population genetics实战

核心模块，

chr.vector，是包含了chr1tochr22的chromosome tag的向量，需要先将其输出到cluster，才能够应用于cluster计算

但是我还存在的疑问：为什么input_df和gene_bed就不export就不会报错了呢？

Note：这里的理解，就好比神经网络模型使用CUDA时，也需要将对应张量加载到GPU上
```
parLapply(cluster_1, chr.vector, df = input_df, fun = function(input_df, chr_number, gene_bed){
  # retrieve output dataset using chromosome tag {1..21} to split the FST dataset
  chr_number_input_df <- input_df %>% filter(chr==chr_number)
  ...
}
```

cluster_1 <- makeCluster(16) # using 16 cores to speed up the FST calculation process.
clusterExport(cluster_1, "chr.vector")  # export chromosome vector the sockets

# 
FST_list <- parLapply(cluster_1, chr.vector, df = input_df, fun = function(input_df, chr_number, gene_bed){
  # retrieve output dataset using chromosome tag {1..21} to split the FST dataset
  chr_number_input_df <- input_df %>% filter(chr==chr_number)
  new_chr_number_input_df = chr_number_input_df %>% left_join(gene_bed, by=c("chr"))
  
  # only save the validated records
  new_chr_number_input_df <- new_chr_number_input_df %>% filter((pos >= start_pos & pos <= end_pos))
  
  # using 0 to replace negative value
  new_chr_number_input_df$weight_FST[new_chr_number_input_df$weight_FST < 0] <- 0

  # calculate FST for each gene
  cal_df <- new_chr_number_input_df %>% 
              group_by(gene) %>%
              summarize(
                Mean_Fst = mean(weight_FST, na.rm = TRUE)
              )
  return(cal_df)
})

# close cluster
stopCluster(cluster_1)